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Abstract 
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1 Introduction 



The problem to include fermionic fluctuations in QCD simulations has been 
in the focus of interest of the lattice community for a long time. An overview 
over the standard approaches to such simulations can be found in |T] . About 
the most obvious idea that one may have, is the following succession of steps: 
First we propose changes of the gauge field by some efficient algorithm that 
fulfills detailed balance with respect to some suitable pure gauge action. 
Then we add a Metropolis accept/reject step. This correction has to filter 
the ensemble such that it becomes governed by the desired action including 
the fermion determinant, whose change enters the accept decision. Since the 
proposal is not guided by the fermions, one may however fear to get sufficient 
acceptance only for tiny rather local changes of the gauge field and to get 
an overall very inefficient algorithm by a succession of such steps. Alterna- 
tive approaches, including the presently most popular hybrid Monte Carlo 
algorithm [21 E] (HMC), therefore use some sort of stochastic representation 
of fermion effects to guide the gauge field evolution at the expense of intro- 
ducing additional noise and having to make many small and expensive steps 
involving inversions of the Dirac operator. 

While fermionic guidance may prove indispensable for many lattice sim- 
ulations, in our opinion there is also some interest to pursue the direct 
Metropolis approach. One reason is the growing interest in Dirac operators, 
where the fermions are coupled to smoothed SU(3)-projected averages of the 
fundamental gauge fields [H El El • Here the dependence of the operators 
on the fundamental fields to be updated becomes so complicated, that their 
change even under infinitesimal moves, which leads to the fermionic force, 
becomes impracticable to evaluate both on paper and CPU. Metropolis on 
the other hand requires nothing but a routine that is able to apply the Dirac 
operator to given fields. A similar situation prevails with Ginsparg Wilson 
fermions, for instance in the form of the overlap formulation p. 

An additional important motivation for the present investigation stems 
for us from our interest in simulations of the Schrodinger functional with 
dynamical fermions fOl El. Here one is also interested in larger /5-values 
implying small physical volume and approximate validity of perturbation 
theory. Then we expect fiuctuations of the determinant to become a small 
(one-loop) effect. Our future hope is to develop a Metropolis algorithm for the 
determinant whose efficiency in the large (3 limit becomes more similar to pure 
gauge simulations based on hybrid overrelaxation. HMC-type algorithms on 
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the other hand remain rather costly also in this limit. If such an algorithm 
can be constructed, it will be interesting to see, if and where there is a cross- 
over in efficiency compared to HMC 

Reliable algorithmic optimizations with dynamical fermions in four di- 
mensional theories are very costly and large scale projects by themselves. 
Therefore — like other researchers — we decided to first undertake a study 
in two-dimensional QED, the Schwinger model ^21 El El El- This simplifi- 
cation allows for clean clinical tests using for instance the precise knowledge 
of the full spectrum of fermionic operators. We are of course aware of the risk, 
that the smoother gauge fields in this superrenormalizable model may teach 
us lessons that do not carry over to QCD. Therefore we plan to soon test the 
algorithms derived here in the four dimensional Schrodinger functional. 

Other efforts to use the determinant directly in Metroplis steps have been 
reported in pTCl for Wilson fermions and in J7] for staggered fermions with 
blocked links. In ^Hl a hierarchical system of acceptance steps has been 
tested. Although interesting, we think it is fair to say that none of these 
projects has led to a strong competitor for HMC for standard actions to 
date. 

Our two dimensional study presented here is organized as follows. In sec- 
tion 2 we set up the notation and our lattice formulation of the Schwinger 
model. In section 3 we investigate the behavior of the global Metropolis al- 
gorithm with determinants evaluated exactly. Our main result is contained 
in section 4 where the stochastic estimation of determinants is introduced 
together with a new class of partially stochastic updates, which is tested in a 
number of applications in section 5. After conclusions two appendices follow 
where we derive exact formulas for the acceptance rate as a function of the 
eigenvalues in an associated generalized eigenvalue problem. The perturba- 
tive solution of this problem is discussed in appendix B. 

2 Model laboratory 

In this section we introduce our formulation of the Schwinger model dis- 
cretized as two dimensional noncompact U(l) gaugefields and Wilson fermions 
Quenched gaugefields are generated by a global heatbath. We work in lattice 
units setting the lattice spacing a = 1. 
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2.1 Formulation of the path integral 

Gauge potentials A^{x) are taken as the primary fields which are integrated 
over all real values. In terms of the field strength 

F^, - A^A, - A,A^, (2.1) 

where is the forward difference, the gauge action reads 

Sg[A]^Wf',. (2.2) 

with /i, z/ = 0, 1. A well-defined path integral on a finite torus of length L 
yields the pure gauge partition function 

Zg^ f DAl[6{j:A,)l[6{A;A,)eM-SG[A]l (2.3) 

II X X 

where A* means the backward difference and DA — H^m^^mC^)- ^~ 
functions fix all modes that do not receive damping by Sg- In addition to 
the gauge degrees of freedom these are two modes corresponding to constant 
shifts of Af^ that we shall come back to. For later use we abbreviate the 
normalized full gauge measure as 

Dpi{A) = ^DA n<^(E^M) ri'^i^;^.) ^m-sg[a]) (2.4) 

G n X X 

and the gauge average as 

{X{A))g = I Di,{A)X{A). (2.5) 

To couple fermions to A^j^^x) in a gauge invariant fashion we choose a 
coupling strength g and form phases 

U^ix) = eMwAi^ix)) (2.6) 

and covariant difference operators 

DJ{x) = U,ix)f{x + il)-f{x) (2.7) 
D;f{x) = f{x)-U;{x-il)f{x-fi). (2.8) 



4 



Now the Wilson operator reads 



Dw^^I{i,{d, + d;)-d;d,} 



(2.9) 



with some choice of hermitian 7 matrices. In terms of 



75 = «7o7i 



(2.10) 



we have the pseudo-hermiticity property 



'Iv = 75-0^^75 . 



(2.11) 



For our algorithmic study we choose periodic boundary conditions for 
the fields that Dw acts on. Had we allowed constant components in 
then we could transform them away by a non-periodic gauge transformation. 
This would however modify the fermion boundary conditions by extra phase 
factors. Our constraint may thus be viewed as a definite set of boundary 
conditions in imposing a finite volume. 

The partition function for Nf flavours of mass m is taken as 



In the following we restrict ourselves to the strictly positive case Nf — 2 
analogous to QCD with only light degenerate flavours. 

Wilson loops constructed from the phases decay in the gauge ensemble 
with an exact area law. The string tension 



is used to eliminate the dimensionful coupling g in favour of the dimensionless 
combination 




(2.12) 




(2.13) 



z — \faL . 

In the pure gauge theory the limit L 
continuum limit at flnite physical volume. 
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2.2 Generation of gauge fields 

Gauge fields distributed with Dfi{A) can be generated by a global heatbath 
procedure or independent sampling. A potential A in the Markov chain is 
followed by A' which, due to the constraints, can be written as 

4 = e^,A>. (2.15) 

The lattice scalar (p is taken as the Fourier transform 

0(^) ~ T i'ip) exp(ip ■ x). (2-16) 
L p 

of independent Gaussian random numbers (j){p) 

0(0) = 0, 4>*{p)=4>{-p) (2.17) 
(r(p)0(g)) = ^5,, (p^O). (2.18) 

Momenta are summed over the appropriate Brillouin zone and depends on 
them periodically. 

It will also be of interest for us to mimic smaller update moves A — » A" 
which are not independent but just fulfill detailed balance with respect to 
Sq- This can be achieved by taking 

A" = cA + sA', c = cos(t7r/2), s = sin(t7r/2), (2.19) 

where the parameter < t < 1 allows to control the step-size. 

In FigC] a few complete spectra of the Wilson Dirac operator are shown 
for several couplings at L = 8 in gauge fields generated according to (|2.4p . 
The high degeneracy of the free spectrum is progressively lifted as z is raised. 
At the same time the spectrum moves away from the origin. 

For each configuration we define an effective critical mass mo to be the 
negative spectral gap 

mo{A) = -mm[Re{spec{Dw))] (2.20) 

such that the smallest real part of all eigenvalues of Dw + mo just reaches 
zero. Its gauge average is taken clS db critical value 

m,{z,L) = {mo{A))G (2.21) 
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Figure 1: Typical spectra of Dw at couplings z = 1, 2, 5, 10 



and by a'^{z,L) we denote the corresponding variance. In Table [T] some 
numerical values are listed. We divide by g"^ as suggested by perturbation 
theory and obtain numbers that vary only slowly with L and z. For the 
spectra in FigQ this implies gaps roughly proportional to in agreement 
with the four particular configurations shown. In our algorithmic study we 
find it appropriate to use these data to choose mass values such that our 
fermions are light and thus dynamically relevant, but heavy enough to not 
suffer from the 'exceptional' unphysical modes known to occur with Wilson 
fermions. 
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L = 8 


L = 12 


L = 16 


z 




rric/g'^ Oclg^ 


fnclg^ Oclg^ 


1 

2 
5 

10 


-0.323(3) 0.088 
-0.314(3) 0.086 
-0.289(2) 0.071 
-0.169(1) 0.020 


-0.382(3) 0.082 
-0.379(3) 0.084 
-0.353(2) 0.068 
-0.262(2) 0.028 


-0.427(3) 0.086 
-0.420(3) 0.081 
-0.391(2) 0.069 
-0.319(1) 0.031 



Table 1: Values for the critical mass and its standard deviation from 1000 
configurations each. 

3 Metropolis with exact determinants 

In this section we use a rather ideal fermion algorithm that is only available 
in our two-dimensional model: global heatbath proposals with respect to the 
gauge action filtered through a Metropolis step based on the exact fermionic 
determinant. With the cost of computing the determinant by standard linear 
algebra means scaling like L^^ this appears prohibitive beyond D = 2. Here 
on the other hand it will prove to be quite feasible up to medium size lattices 
and will provide a rigorous upper bound for the acceptance rates achievable 
with stochastic techniques. 

3.1 Exact acceptance rate 

In equilibrium for the ensemble ()2.12j) the acceptance rate q for proposals 
with the pure gauge global heatbath described in the previous section is 
given by 

(3.1) 

where A and A' enter into Dt^ and D'^^. In a more symmetric form this 
reads 

q = ^ I Dn{A) J Dfi{A') mm {\det{Dw + m)\^,\ det{D\y + m)\^) 



!Dfi{A)JDfi{A')\ 


det{Dw + TTi) 


2 _ 


det{D'^ + m)\'^\ 


JDf,{A)JDfi{A'){ 


det{Dw + m) 




detiD'^ + m)\^) 



and q clearly obeys < g < 1. Any nontrivial dependence of the determinant 
on the gaugefield reduces the acceptance. 
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To estimate q we generate a large number of independent gaugefields 
with the pure gauge measure and compute the fermion determinant for each 
of them. Let di,i = 1, . . . , iV be the resulting successive values of | det {Dw + 
m)p). Then q may be estimated by 



(3.3) 



TV 



For larger A^- values the double sum is best evaluated by using a sorting algo- 
rithm (with cost only oc A^ln(A^), see jTHj, for instance provided in Matlab), 

{d<} = sort({rf,}), (3.4) 

where the sequence {df} consists of the same members as {di} but reordered 
such that df < df < . . . < d"^. Now the acceptance is written as 

- - (3.5) 



with the weights 



2(N-i) , , 

For the error estimation one of course has to take into account that the 
O(iV^) terms in the numerator of (|3.3|) that are effectively summed by (j3.5|) 
are not independent. 

It turns out to be very successful to make a Gaussian model for the 
distribution of the fermionic action in the gauge ensemble 

v{E)=<b{E-SF)>G. = -21og(|det(Dw' + m)|), (3.7) 

by setting 

Once this Ansatz has been made, its free parameters and, more im- 
portant, S can also be estimated numerically from the observed values di by 
determining mean and variance of Sp. Within the model the acceptance rate 
then follows, 

_ ^ /+- dE v[E) exp(-E) /f^ dE' v[E') 
' \l^dEv{E)^^v{-E) ■ ^'-'^ 
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With the expression (|3.8p for v one sees that q is independent of £"0 as it 
should (irrelevance of a constant in Sp). It is convenient to extract q from 
the distribution z/(A) for the energy difference I\ = E' — E 

\ r+oo r+OD 

i){A) = - dEu{E)exp{-E) dE' u{E') 6{A - E' + E) (3.10) 

Zj J — CO J 00 

with 

/+00 
dE u{E)exp{-E), (3.11) 
-00 



for which we obtain 



In terms of u we evaluate 

/+00 2 
dA z>(A) min(l, exp(— A)) = —= / duexp{—u^) = erfc(S/2). 
-00 yvr JE/2 

(3.13) 

We performed a series of (quenched) simulations on L = 8, 12, 16 lattices 
with z = 1,2,5, 10 (compare Table where we determined and stored the 
complete spectra of D]y. These data may be used to construct fermion 
determinants and acceptance rates for a whole range of masses and in this 
fashion we produce Fig|21 It demonstrates that in our model, at least with 
the exact fermion determinant, a simulation based on global Metropolis steps 
is feasible over a wide range of masses. Note that this even refers here to 
maximally large quenched gauge move proposals. Moreover we verify the 
validity of the Gaussian model (circles). 



3.2 Enhanced acceptance rate 

So far we have considered the generation of potentials with the pure gauge 
action So and the incorporation of the determinant in an accept/reject step. 
The same ensemble may also be produced by a different split of the total 
action and it is conceivable that this leads to an enhanced acceptance rate. 
It corresponds to a simple version of ultraviolet filtering fHl 1201 1^ by mod- 
ifying the fermion action just by the plaquette term. Apart from the gauge 
action contained in the measure Dfi{A) we include another component 

Z= [ Dij{A)exp[{l-a^)SG]\det{Dw + m)f (3.14) 



10 



0.8 
0.6 

cr 

0.4 
0.2 




z=1 



z=10 



C C 



6 , 7,, 8 9 10 
(m-m )/a 



Figure 2: Acceptance rates versus mass for L = 16. Numerical values from 
are shown as crosses and circles correspond to the model using (jH.lHj) 
with the observed E. 



which in simulations is combined with the determinant in the Metropolis 
test. By rescaling one easily sees that this ensemble is equivalent to the 
standard form of the gauge action with an effective coupling of strength g/a 
in ()2.6|) . Whenever the highest acceptance, for fixed g/a, is reached at a 7^ 1 
the extra term has enhanced the acceptance and is hence advantageous. 

In simulations we perform the expensive evaluation of the determinant 
with one fixed value of g in (j2.6|) and then compute the acceptance for many 
values a in the additional term. In this way we construct lines in a graph of 
q versus z, where z refers to the effective coupling. With a number of such 
lines, we shall see which one gives the highest acceptance for a z at which we 
wish to simulate. 

The Gaussian model generalizes by assuming a joint Gaussian distribution 
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for Sg, Sp with a 2 x 2 covariance matrix given by connected correlations 

<SgSg>'g = VGG = \{L'-l) (3.15) 
< SpSp >G = vff (3.16) 
<SgSf>g = vgf = vfg, (3.17) 

of which the first one is trivial due to the Gaussian gauge action. 

Instead of Sp alone, acceptance is now controlled by the combination 
E = Sp + (a^ — 1)5*^ whose variance S^, given by 

S2 = _ 1)2^^^ + + 2{a^ - 1)vgf (3.18) 

may now be used in ()3.13|) to estimate q. Within the model, we can easily 
vary a continuously. 

In FigI21 we see a number of approximate acceptance trajectories con- 
structed in this way for L = 16, m = —0.1. For each of them 1000 gauge 
configurations were used and exponentiated in the determinants with cou- 
pling z corresponding to the crosses. The parameter a within the Gaussian 
model was then varied in the interval [0.8, 1.2] producing the acceptance tra- 
jectories as functions of z = z/a. We clearly see that ultraviolet filtering pays 
off. At weak coupling the acceptance is high without it, and correspondingly 
less can be gained. For the case with the cross at z = 10 we now confront the 
Gaussian predictions for the a-dependence of q with numerical values as in 
fl3.3|) where now di = exp[— S'i? — (a^ — 1)Sg] enters. This is shown in Fig 01 
for an ensemble of 1000 gauge potentials and confirms the model also in this 
more general setting. 

4 Partially stochastic global acceptance steps 

Metropolis steps based on the availability of exact fermion determinant ratios 
are of conceptual interest but hardly lead to efficient algorithms on large 
lattices and in four dimensions. An approximate stochastic estimation that 
nevertheless maintains exact detailed balance can be defined and at first 
seems more promising [TTl 122 . We found however that this can easily lead 
to negligibly small acceptance at smaller masses. Therefore we developed a 
more general Metropolis correction step based on an only partially stochastic 
estimation of determinant changes (PSD) which is complemented by a few 
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Figure 3: Enhanced acceptance rate in the Gaussian model for L = 16,m = 
—0.1. The crosses correspond to a = 1. Each of them extends to a curve 
according to (j3.13j) . (j3.18|l . as a is varied in the interval [0.8, 1.2]. 



exact eigenvalues [see [21] for related ideas in connection with reweighting 
corrections to the polynomial hybrid Monte Carlo algorithm] . In this section 
we shall thus derive a family of (exact) algorithms which contains those based 
on exact and on fully stochastic determinants as extremal cases between 
which we interpolate. 



4.1 Stochastic estimation of fermion determinant ra- 
tios 

The stochastic estimation of the determinant reduces the problem to solving 
linear equations with the fermion matrix at the expense of statistical errors 
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Figure 4: Acceptance rates for L = 16,m = —0.1. The curve corresponds to 
the Gaussian model based on the simulation z = 10, other data with errors 
are direct Monte Carlo checks. The circle is at a = 0.950 that minimizes 

dsmi). 



on top of the main Monte Carlo process. The fundamental formula is 

D[r]]p{Mr]) = \ det(M)rl (4.1) 

Here we integrate over a complex valued spinor field with the measure 

dlie{ri)dlm.{r]) 



TT 



(4.2) 



for each site x and spinor component a. For the normalized probability 
distribution p, 

' D[7]]p{7]) = l (4.3) 
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we often take a Gaussian^ 



p{ri) = exp —rj'^T] , (4.4) 

but keep the formulas more general where we can. The determinant appears 
as a Jacobian in ()4.H) for arbitrary p. We may also write an unbiased estimate 
of the determinant as 

with an average < . >^ over the random field rj only. 

For an acceptance step from a given field A to a proposed new A' with 
associated operators D^, D'^y we now form the "ratio operator" 

M = {D'y^, + m)-\Dw + m), (4.6) 

and in terms of this matrix we could stochastically accept with the probabil- 
ity 

wo{A, A') = min[l, p{Mri) / p{ri)] , (4.7) 

where the dependence on r] and the choice of p is left implicit. For the reverse 
transition, A A', we find M ^ M~^. Therefore 

iMAA'))^ ^ JD[7^]mm[pir^),piMv)] ^ , . .^.,-2 o^ 
{wo{A',A))^ JD[r^]mm[p{r^),p{M-^v)] \ ^ '\ ^ ' ^ 

shows detailed balance. The last equality follows by changing variables 
T) Mr] in the integral in the numerator. These steps constitute the fully 
stochastic algorithm that we are going to generalize to PSD^. 

The acceptance rate is always smaller than (|3.1|) due to the inequality 

D[r]]p{r]), f D[r]]p{Mr]) 



D[ri]mm[p{ri), p{Mr])] < min 

= min(l, |det(M)|-2) (4.9) 

Hence, the exact acceptance rate in this context is something like the ideal 
"Carnot" efficiency which we cannot reach but which we may also not want 



^The sum over x, a is included in the scalar product 77^77 here and in the following. 
^It would most probably be advantageous to use a preconditioned operator here. In 
this study of principles we however avoid this complication. 
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to miss by too much. The question may arise, if the acceptance rate may 
be increased by averaging over several random r] fields^. If we average under 
the min function, A^^ = 1 seems to be the only finite value for which detailed 
balance can be shown. Averaging outside of min seems correct but would 
not raise the average acceptance. 

An expression for the stochastic acceptance rate with distribution p is 
given in analogy to (j3.ip by the integral 

qo = ^ I Dfi{A)\det{Dw + m)\^ I Dp{A') {wo{A,A')), (4.10) 
= ^1 Dp{A)\detiDw + m)\^ I Dp{A') J D[r]]mm[p{r]), p{Mr])] 

= Dfi{A) I Dp{A') J D[r]]mm[p{{Dw + mr'v),PiiD'^ + m)-'v)] 

A naive Monte Carlo estimation of the last expression for go with p = 
exp{—rj^ri) is not practical due to the very strong fluctuations of the inte- 
grand. We shall however be able to perform the ?]-integrations exactly in 
this case, which may be viewed as the construction of an improved estimator 
(same mean value, smaller variance) for go in terms of generalized eigenvalues. 
We work out the dependence of 

{wo)^ = J D[r]]mm[p{r]),p{Mrj)] (4.11) 

on the spectrum {Aj} of M^M with i = 1, . . . ,n = 2L^ for the choice p = 
exp(— T]''"?]). Performing the above integration in the basis of orthonormal 
eigenvectors of M'^'M with components Zi we find 



K).=n(/ 



dKeZidlmZj ' 



mm 



Changing to polar variables in all the complex planes we get 



{Wo)rj = n 



mm 



Ui 



(4.12) 
(4.13) 



In appendix |X1 this integral is evaluated exactly yielding 



i jjLi - 



^This would be possible, if the determinant were incorporated by a reweighting instead 
of an acceptance step ■ 
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where the special case of S being the empty set suffices here (compare ()A.1|) ). 

It is clear from (j4.13p that eigenvalues Aj = 1 are irrelevant for the ac- 
ceptance. Approaching this limit for one of them in ()4.14|) one indeed finds 
it to drop out and one is left with the formula for n — 1 eigenvalues. If we 
now consider as an example the case of A2 . . . A„_i differing from one only 
negligibly and only one remaining pair with A„ ^ 1 > Ai then we find a 
small acceptance 

K).^^. (4.15) 

It remains small even for AiA„ = 1 = det (M^M), when we have 100% non- 
stochastic acceptance. This simple example demonstrates how the stochastic 
acceptance rate degrades if a determinant ratio of order unity arises from 
compensations between the eigenvalues of the squared ratio operator M'^M. 
We are hence motivated to have a closer look at such spectra in our model. 

4.2 Spectrum of random quotients of Dirac operators 

In practice we find the spectrum {Aj} of M'^M by solving the generalized 
eigenvalue problem 

{Dw + m){Dw + m^x = KD'w + m){D'y^ + m)^x ■ (4.16) 

A general observation about the spectrum of M^M, at least for not too strong 
coupling, is that there are two low and two high almost degenerate eigen- 
values separated from the remaining ones. In the bulk of the spectrum the 
eigenvalues are close to one which is easy to understand, since all eigenvalues 
would be exactly one if the two gaugefields entering into M would be equal. 

A few experiments show that the separation of the extremal eigenvalues 
becomes more pronounced as the mass is lowered. In the same limit and at 
small coupling they are of the order g'^ and g^'^ respectively. This is illus- 
trated in Fig|31by the example of the quenched spectra on L = 8 for several 
combinations m, z. They are labelled by A-E with the corresponding values 
contained in TabI21 in Sect. 15.11 The plot shows the complete generalized 
spectra on a logarithmic scale. For the spectra labelled by A-D one clearly 
sees that the bulk is between 1/2 and 2 with two eigenvalues close to g'^, g~^ 
on either side, while for the more strongly coupled case E this segregation is 
lost. 
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Figure 5: Typical complete generalized spectra {Aj} for values m, z found in 
TabEl Vertical bars show the scales and g~^. 

To understand how such a spectrum arises, we consider two hermitian 
operators H and H' and the generalized Ritz functional (or Rayleigh quo- 
tient) 

The vectors '0* for which /i(V') is stationary satisfy the generalized eigen- 
value equation Hip^ = XH'ip^ with A = /i(^/'*). In particular the smallest 
generalized eigenvalue Amin has the property 

/i(^) > Amin for all?/'. (4.18) 

In our case H = {Dw + fn){Dw + m)^ and H' = {D\^ + m){D'yy + m)^ The 
typical situation is encountered if we set m = and assume that H and H' 
are positive for nonvanishing coupling (note that rric < 0). 
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Instead of considering directly the generalized eigenvalue problem (j4.16|) 
we use the relation H = DwD\^ = {Dw'y5y and first study the spectrum 

-Dh^75V^ = k^ (4.19) 

in perturbation theory. By expanding t/^ = 1 + igA^ — + 0{g^) we 

obtain an expansion of the Wilson Dirac operator 

= + gD^^^ + g'D^S + 0{g')) . (4.20) 

We are interested in the perturbation of the eigenvalue k^^^ = of the free 
operator D^)'-f5. It is doubly degenerate with spatially constant free eigen- 
functions We set up expansions 

i^a = AL'^ + 0((7'), (4.21) 

i;^ = ^(^)+gi;W+g'i;(^) + 0{g'), « = 1, 2 . (4.22) 

The perturbation D^) is proportional Afj,{x) which has no zero momentum 
component. Hence matrix elements of this operator vanish between states of 
equal momentum and in particular between ijj^^ . This implies the vanishing 
of K^^\ Suitable zeroth order eigenfunctions ijj^^ are linear combinations of 
ijj^'^ which we determine later. The first order equation for the eigenvectors 
is 

= « = 1,2. (4.23) 
In second order the 2x2 matrix 



# [D'w' - D\;^>D^^>-D\^>\ 75^r , kj = ± (4.24) 

has to be diagonalized. Here the eigenvectors ip^^ get determined together 
with eigenvalues n^'^ which lift the degeneracy. 

Now we consider the Ritz functional /U (■?/') (j4.17|) setting ip equal to ipa = 
i^a'^ + g^a^ + 0(5'^)- Using the relation (OT^ we get 

i^lDwDl^ij^ = (42))2^4 ^ Q(^5) ^ (4 25) 

i^iD'^D'U. = ^^'liD^w^-D'^^HM^^g' + Oig'). (4.26) 

From ()4.18|) we derive an upper bound for the minimal generalized eigenvalue 

Amin of dnni) 

Amin < const X g^ (m = 0) . (4.27) 



i(2) n(i)n(0)-in(i) 



(0) 
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By considering the inverse Ritz functional fi~^{ip) and setting ip equal to the 
eigenfunctions ip'a of -^vkTs can derive in full analogy to the above a lower 
bound for the maximal generalized eigenvalue Amax of ()4.16p 

Amax > const' X g'^ (m = 0) . (4.28) 

A detailed analysis of the perturbation expansion of the generalized eigen- 
values themselves, which confirms the variational arguments given here, is 
deferred to appendix iBl 

4.3 Partially stochastic estimation of determinant ra- 
tios 

We saw that a few extremal eigenvalues in M^M can ruin the stochastic ac- 
ceptance rate even if the relevant determinant ratio is close to unity. More- 
over, at least in our model, these kind of spectra are really common. We now 
develop a mixed strategy treating the bulk of the eigenvalues stochastically 
and a few special ones exactly. Of course, detailed balance has now to be 
demonstrated for the combination. For the time being we content ourselves 
with proofs for p{ri) = exp(— r^"''?]) only. More general cases may well be 
possible. 

We assume, that for each spectrum we can identify a set S such that 
{Xi\i G S} consists of the s = jS"! extremal'* eigenvalues of M^M (the same 
number of large and small ones), which we treat deterministically. The as- 
sociated eigenvectors {4>i\i G S}, that we choose to be orthonormal, span an 
s dimensional subspace that we characterize by a projection operator 

P = ^0,0|^P(AA') (4.29) 

Here we remind ourselves that P, via the Dirac operators Dw, D'w, depends 
on a pair of gauge potentials A, A'. Due to the symmetric inclusion of large 
and small Aj a spectral analysis of the inverse operator (M'''M)~^ with in- 
verted eigenvalues would lead to the same projector (for the same M). 

Anticipating a discussion of detailed balance we are interested in informa- 
tion about the relation between P{A, A') and P{A', A). As mentioned before, 
going to the reverse process [A ^ A'), M changes to M~^. Hence in this case 

^We assume this construction to be unique, i. e. no degeneracy at the boundaries of S. 
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we are concerned with the eigenvalue problem of M~^^M~^ = {MM^)~^. Its 
eigenvalues are reciprocal to those of M^M but the eigenvectors and hence 
the associated projectors are different, P{A',A) ^ P{^A^A!\ The problems 
are however related by a unitary transformation 

MM^ = U^M^MU, U^U = 1 = f/f/t . (4.30) 

Explicitly, U may be written as 

U = M^MM^)-^/^ . (4.31) 

The same unitary transformation relates the eigenvectors of the two problems 
such that 

P{A\ A) = U^P{A, A')U . (4.32) 

holds. 

Returning to the forward problem we can factorize 

M^M = (P + PM^M) (P + PM^M) (4.33) 

where we introduced the complementary projector 

P = 1 - P (4.34) 

and used standard properties of projectors and commutativity PM^M = 
M^MP, which is easily seen in the spectral representation of M^M. While 
both factors in ()4.33|) are nonsingular operators in the full domain, they have 
a block structure with unit operators in the subspaces of P and P. The exact 
"small" determinant of the first factor is 

det(P + PM^M) = n Ai. (4.35) 

i&S 

A stochastic estimator for an acceptance step with the second factor alone 
would be given by (see ()4.7p ) 



mm 



1, exp(— r/' 



As the true partially stochastic acceptance criterion Ws{A,A') to accept a 
proposed A' given an 'old' A we now propose^ 



WsiAA') 



mm 



1, n exp{-7]^P{M^M - 1)Pt]) 



(4.36) 



^One could also think of two separate successive Metropolis steps, but this leads to 
smaller overall acceptance rates. 
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To prove detailed balance we start from 



1, n exp{-ri^U^PU{{MM^)-^ - Ip^PUrj) 
ies 



(4.37) 



Ws{A' , A) = min 
In 

{ws{A',A)),= (4.38) 

Dlr]] min exp{-r]^r]), Y[ Xiexp{-r]P\P + P{M^ M)-^)Ur]) 
I ieS 

a change of variables 

= [/t(p + PM^Mf'^ri' (4.39) 



yields the desired result 

{ws{A,A')), 



{ws{A\A)) 



n A,^Met(P + PM^M)-i = |det(M)|-2. (4.40) 



With the above formalism we have an algorithm which avoids the low 
stochastic acceptance from extremal eigenvalues if their product is of order 
unity. It remains to discuss how to compute the required eigenvalues and 
-vectors. We do not attempt a detailed discussion here. An obvious idea is 
however to generalize the method used by the ALPHA collaboration in the 
past to obtain low and high lying eigenvectors of the ordinary problem. It 
is based on minimizing a Ritz functional by a conjugate gradient technique 
|26j . The relevant functional for the generalized problem is 

MX) = li^r't'tt 

\{D^ + m)xV 

This functional is extremal at generalized eigenvectors and its value there 
is one of the generalized eigenvalues which coincide with those of M'^M. 
Different generalized eigenvectors Xi ^^e not orthogonal, but one may show 
that x'i^Dw + m){Dw + 'mYxj = if Aj 7^ \j. In fact, {Dw + myxi ^^re the 
(unnormalized) eigenvectors of M^M and {D'^r + fnj'^Xi ^i^re those of MM^ 
After finding the absolute minimum of /i at xi o^ig i^^iy then search for 
the minimum in the space orthogonal to {Dw + Tn){Dw + ^)^Xi to find X2 
belonging to the second smallest eigenvalue A2. For the largest eigenvalues 



22 



one may proceed analogously or just exchange numerator and denominator 
of /i. Note, that we do not need a solver during the minimization, that is 
no intolerable nested iterations. A solver is needed however in the stochastic 
part to apply M or M^M to Prj. Due to the filtering of the random vector 
through the projector this should however be a well-conditioned problem 
and take few iterations. An alternative method to construct the required 
eigenvalues and vectors could be Lanczos techniques j2Zl which are known to 
first converge for the extremal eigenvectors needed here. 

5 Numerical experiments with PSD 

5.1 Acceptance rate for PSD 

We now consider the partially stochastic acceptance rate 



The 'observable' {ws)n in this double pathintegral depends on the generalized 
eigenvalues Aj and on their division into the deterministic subset 5* and the 
stochastic one S. As limiting cases it includes the fully deterministic and 
fully stochastic evaluation. For a numerical estimation of Qs we produce 
pairs of configurations A^, A'^ and determine for all of them, at some value of 
z, dk = \ det{Dw + m)p, generalized eigenvalues Aj and from them estimates 
flfc = {ws)r] according to (jA.HIA.lGll . Then we have the Monte Carlo estimate 



In Figini it is shown how the average acceptance rises if starting from the 
fully stochastic case s = more and more eigenvalues are included in S. 
Parameters are L = 8,2; = l,m = 0.0125 and 1000 pairs of configurations 
are sampled. Some further examples can be found in Tabj2I We note that 
there is a roundoff problem in the straightforward evaluation of ()A.16j) due 
to cancellations and significance loss. With standard double precision accu- 
racy some clever recombination of terms would be required before using the 
formula much beyond L = 8 with n = 128 eigenvalues Aj. 

At this point we also study the dependence of the acceptances on the 
'size' of the proposed moves. Global heatbath steps with respect to the 
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Figure 6: The acceptance rate vs. the number s of nonstochastic generahzed 
eigenvalues (parameters B in Tab|2I). 

gauge action that we have used so far are of maximal size. Lowering from 
one the parameter t that we have introduced in ()2.19p allows us to mimic 
smaller steps. In Figl7| we see the expected dependence of the acceptances. 
In particular we see that partially stochastic acceptance steps allow for much 
larger moves without excessive rejections. 

Also in the partially stochastic case we have found a successful Gaussian 
model for the acceptance rate. Our starting point here is ()5.1j) which we 
break into two steps 

z>,(A) = (( ( 5[A - 5: In A, - ri^PiM^M - l)r]] ), )) (5.3) 

ies 

with 

{{O)) = i| Df^{A)\det{Dw + m)\' I Dix{A')0 (5.4) 
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L z 


m 


1 


% 




label 


8 1 


0.0250 


0.837(7) 


0.0141(2) 


0.460(4) 


A 


8 1 


0.0125 


0.734(11) 


0.0029(1) 


0.425(6) 


B 


8 1 


0.0050 


0.602(15) 


0.00061(2) 


0.368(8) 


C 


8 2 


0.0125 


0.634(14) 


0.0020(1) 


0.130(3) 


D 


8 4 


0.035 


0.819(8) 


0.00084(4) 


0.0083(2) 


E 



Table 2: Values for acceptance rates for the deterministic (g), stochastic (go) 
and mixed case (^4) with four extremal eigenvalues in S. 

and 

/+00 
c/Ai>s(A) min(l,exp(-A)). (5.5) 
-00 

As a model we assume t's(A) to be a Gaussian of mean and width bg- For 
it Qs evaluates to 

Qs = - erfc I I H — erfc exp(6^/2 — mo) . (5.6) 

In the expressions for the mean and width of the A distribution the rj inte- 
grations can be performed and we find 

= ((^ln\ + E(A.-l))) (5.7) 

In TabOthe last column is replaced by (0.452,0.419,0.364,0.140) in the Gaus- 
sian model for A-D, while the agreement is worse for case E with its different 
pattern of fluctuating eigenvalues. 

5.2 The relevance of gaugefixing 

We now turn to the question of gauge fixing. Up to now all gaugefields were 
generated in one and the same completely fixed gauge according to ()2.4|) . 
Stochastic acceptance steps depend on the operator M of ()4.6|) . If only 
p{ri) is invariant under x-dependent phase transformations of the random 
field ri{x) then all acceptances are invariant under a change of gauge, that 
is the same transformation applied to A and A'. If however the update 
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Figure 7: Acceptances rates versus stepsize t for L = 8, m = 0.0125, z = 2. 

proposal includes a more or less random gauge move of A' relative to A, this 
is equivalent to only gauge transforming A'. In this case the nonstochastic 
acceptance probability involving the ratio of gauge invariant determinants 
as in ()3.H) is unchanged. This is however not true for the eigenvalues Aj of 
M^M and of the generalized eigenvalue problem ()4.1(jj) . We looked at the 
generalized spectrum for pairs Afj,,A'^ + Af^j for fixed A, A' and a number of 
randomly chosen gauge functions 7(0;). For nonvanishing 7 the eigenvalues 
Aj show much more variation and smaller values of -F(Aj) in ()A.1|) result while 
the ratio of determinants Hi remains unchanged as it has to. 
A non gauge-fixed simulation is achieved by replacing in (j2.4p 

T[S{A;A,) exp (-^E(A;^m)') (5.9) 

such that the value ^ = corresponds to the previous case. It is not difficult 
to show that a global heatbath step for this measure is effected by a gauge- 
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fixed step as before followed by a gauge transformation with a gauge function 
7 = ^0. Here a new independent random field is produced obeying (j2.171 
I2.18p . Acceptance rates as functions of the gauge parameter are shown in 
FiglSl We notice the dramatic decay of the already small purely stochastic 
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Figure 8: Acceptance rates versus gauge parameter ^ for L = 8, m = 
0.0125, z = 1. 

acceptance. 



5.3 Some unquenched simulations with PSD 

As a practical application we performed a few simulations on lattices L = 24 
with z = 1 and m = 0. By looking at TabQwe observe that for this value 
of z the value of ac/g"^ seems to depend only weakly on L whereas —rric/g'^ 
increases monotonically with L. We have good reasons to expect that at 
these parameters no 'exceptional' fields with zero eigenvalues of the Wilson 
operator D\y + m are proposed. 
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We compare PSD simulations at s = 4 with exact determinants (equiva- 
lent to s = 2L^). At s = 4 we employ both global heatbath {t = 1) proposals 
and smaller update moves with the parameter t in ()2.19|) set to 1/2. We 
measure the 'pion'-susceptibility in the unquenched ensemble 



where Tr refers to both Dirac indices and space. Before integrating out the 
Grassmann valued fermion fields x can also be interpreted as a susceptibility 
of the density ilji{x)'j5ipj{^) for two flavors ipiji = 1,2. We take x h^re as a 
simple correlation function with significant contributions at the scale of the 
meson correlation length. 



Table 3: Results from simulations of 4900 measurements each on L = 24 
lattices with z = 1 and m = 0. 

In TabOlwe list our results for the acceptance rates and the susceptibility 
X together with its integrated autocorrelation time Tmt.x; which is in units of 
(global) gaugefield updates. The stepsizes t = 1/2 and t = 1 are equally 
efficient here, since the smaller but equally expensive steps are precisely 
balanced by the larger acceptance. An optimization is not attempted in 
this model study. It will look different, if the stepsize is also tuned by the 
number of modified links and not just by the amount of change. 

In the simulations with PSD we also looked at the distribution of the 
action difference A that appears in ()5.3|) . In FigElwe show histograms with 
the distributions of the two components contributing in ()5.3p for the two 
choices t = 1 and t = 0.5. From the measured values of A we extract its 
mean rris and variance b^. For t = 1 we get = 2.82 and = 5.90, for 
t = 0.5 we get = 0.92 and b^ = 1.88 both a factor three smaller than 
for t = 1. Inserting these values in ()5.6|) we obtain acceptance rates in the 
Gaussian model of 24% for t = 1 and of 50% for t = 0.5, in perfect agreement 
with the acceptance rates directly observed in the simulations. 




(5.10) 



s t qs X Tint,x 



2L2 1 0.61 983(12) 2.4(3) 
4 1 0.24 1005(15) 4.4(8) 
4 1/2 0.50 997(14) 3.8(5) 
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Figure 9: Distribution of J2ies^^^i (left) and ini'^P{M'^M — 1)7] (right) con- 
tributing to u in (|5.3p . 

6 Conclusions 

We found that two-dimensional QED can be simulated over a wide range 
of parameters by a combination of pure gauge update proposals combined 
with a global Metropolis step including the exactly evaluated determinant 
for two flavours of Wilson fermions. Employing a coupling for the proposals 
that differs from the physical one can raise the acceptance rate further and is 
an example of ultraviolet filtering. This algorithm provides an upper bound 
for the acceptance rate of a class of stochastic techniques. We studied such 
methods where we replace the ideal exact determinant ratios by less expensive 
stochastic estimates whose main cost are inversions of the Dirac operator. 

It turned out to be both necessary and possible to find a mixed form of 
stochastic estimates combined with the exact incorporation of a few critical 
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modes which otherwise spoil the acceptance probabihty. These modes are 
given as extremal eigenvalues of the generalized problem defined by the Dirac 
operators in the present and the newly proposed gaugefield. This was coined 
PSD algorithm as it is based on the partially stochastic determinant. Gauge 
fixing proved to be advantageous in this case. If omitted, the update propos- 
als include random relative moves along the gauge orbit which are punished 
with enhanced rejection. For the physically smaller volumes we found the 
acceptance rates to be well described by assuming Gaussian distributions for 
those parts of the total action that enter the Metropolis decision. The rates 
are then given by error functions of the observed mean and variance. With 
these insights we plan to investigate similar update schemes for four dimen- 
sional QGD, especially at small and intermediate physical volume, and hope 
to report on this in the future. 

Acknowledgements We wish to express our thanks to Burkhard Bunk, 
Anna Hasenfratz and Rainer Sommer for helpful discussions, and e-mail ex- 
change with Tony Kennedy is acknowledged. All calculations were carried 
out in Matlab employing its built-in linear algebra tools. 

A Acceptance rate for given spectrum 

In this appendix we evaluate the integral 

= n ( / "^i^ exp(- n exp(- ^ XiUi) . (A.l) 

ies ^"^^ ^ ies ies ies 

Eigenvalues < Aj < oo,i = 1 . . .n enter here and S is a subset of their 
indices which may also be empty. The set S comprises the remaining indices. 
The part, where the first argument of min is minimal, is given by 

fi>^^) = I[(r du^) exp{-J2m) 4E(1 - - (A-2) 
ieS ^ ° ^ ieS ieS 

where we abbreviated 

C = ^lnA,. (A.3) 
ies 

By rescaling integrals the full result is obtained as 

F{Xi)=f{X^) + {llK')fiK'), (A.4) 
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a form which immediately reflects detailed balance in the form 

To evaluate / we use the Fourier representation of the step function 

^ (A.6) 

where u is real, e positive and infinitesimal, and the integration runs along 
the imaginary axis. Proof: residue theorem. With this inserted into ()A.2|) 
the M-integrations factorize and can be carried out to yield 

/+ioO rl-/ p — zC" 1 
1^ V n 1 — rr^v (A-^) 
-ioo 2111 z + e J-^^ 1 — z[l — Xj) 

Assuming non degenerate Aj also this integral is evaluated by the residue 
theorem. The contour can be closed with negligible contribution in the right 
or in the left half of the complex plane depending on whether C is positive 
or negative. Performing the integration for both cases we find 



J:^esx<l e^-^ n,e5,,v. jE^ for C> 
/(A.) = <! ^ (Ai 

1 - E.e5,A.>i e^-^ n,,5,,v. fe^ for C < 



In the special case C = either closure of the contours is legitimate and 
the two expressions coincide. The implied identity, which holds for arbitrary 
non- coinciding Aj, 

n n 1 

i=l i¥'j=l * 

can in fact also be verified purely algebraically. To this end we recall the 
Lagrange polynomials, which are familiar from interpolation and numerical 
integration, 

- 

with the property 

hiX,) = S,,. (A.ll) 
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We may also introduce a matrix A with elements Ojj by expanding 

ki\)=J2a,,\^~\ (A.12) 
j 

and in particular its last column is given by 

j^i - ^3 

The famous Vandermonde matrix V has matrix elements 

v,k = ihy-' (A.14) 

and in terms of it ()A.11|) says that A and V are inverse to each other. A 
particular consequence is the relation 



which also implies ()A.9|) upon binomial expansion of the numerator. 

It now remains to combine the two contributions to F in ()A.4|) . A short 
calculation leads to the result 

+ E,,e5,A,<i(l - 1/A.) n,g5j^i ^ for C> 

1 - E.65,A.>i(l - 1/A.) n,,5,,y. ^ for C < 

(A.16) 

Note that here the exponentials are both damping and may also be written 
as ^ ^ 

e— = n(Afc)^. (A.17) 

Zees' 

For S* = we recognize the deterministic acceptance min[l, Hfc A^^]. For the 
fully stochastic case 5 = we may use identity ()A.9|) and its companion 



E l(A. - 1)"-^ n = n (A-18) 



to obtain (jHHj). 
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B Perturbation theory for generalized eigen- 
values 



B.l General expressions 

We want to solve the generalized eigenvalue problem 

Hx = XH'x (B.l) 
by finding A as stationary values of 

The operators H = , H' = h'^ (all hermitian) possess expansions 

h = /lo + E/^fc (B-3) 

k>l 

h' = ho + ^g% (B.4) 

fc>i 

and similarly for H, H'. We are first interested in the case that ho has a zero 
mode 

/io0o = 0, 0^00 = 1, (B.5) 
which in perturbation theory is only lifted in second order because we assume 

(f)lhi(f)o = (f)lh[(f)o = 0. (B.6) 

In this case there is an eigenvector = 0o + ^701 + 0{g'^) of h, 

h(j)^K(j), K = /«2/ + 0(/). (B.7) 

with 

/io0i + hi(f)o = (B.8) 

and 

K2 ^ 4{h2 - h{ho)-^h)(f)o. (B.9) 

We next want to show that to leading order in g the vector is also a 
stationary point of r with the corresponding generalized eigenvalue A = 
r(0)(l + O(^)). 
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We digress here to a simpler case to show the structure of the argument. 
Assume we want to know an extremal value of r{x) = p{x)/q{x) close to an 
already known extremum p'(a;*) = with a small value of p{x^)/q{x^) = e. 
Then we set x = x^ + y, expand, and find that the extremum of r is at 

y = e^^+0{e') (B.IO) 

with value r(x* + ?/) = e + O(e^) unchanged to leading order of the expansion. 
Returning to the real problem we set 

X = ^ + R, i?V = 0. (B.ll) 

Stationarity of r implies 

Hxii\x) = H'xii{x) (B.12) 
It is not too difficult to see that with our Ansatz we get 

Mx) = g'4 + 0{g') + \hR\^ (B.13) 
h'x = giK~h)<Po + 0{g^) + h'R (B.14) 

There is a solution to the above equations of the form 

R = g^Y.9'Rk. (B.15) 

fc>0 

This implies = + 0{g^), = + 0{g^) and 

A = X,g^ + 0(g^) = r{<l)) + Oig') (B.16) 

- /ii)0or 

Using ()B.12|) the leading term Rq is determined by 

hoRo = (K - ^i)0o A2, 0i)i?o = (B.18) 

with the last condition following from ()B.11|) . 

There is another generalized eigenvalue associated with the smallest eigen- 
value of h'. Since the problem has the exact symmetry H ^ H', A ^ 1/A it 
is clear that the corresponding eigenvalue is A' = A'_2fi'^^ + 0{g^^) with 

,_^JSK^ (B.«) 
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where n'2 is as 1^2 in 

(PT|) but now referring to the lowest eigenvalue of 
h' . The respective shift R' of the eigenvector 0' has a leading component 

We now come to the regular generalized eigenvalues associated with the 
non-nuUspace of /io, i-e. with zeroth order different from We expand 

X = Y.a'Xk, hoXo^O (B.20) 

fe>0 

A = i + Y.a'h (B.21) 

k>l 

and insert this into (IB. II) to derive for the first two nontrivial orders 



A2 = ^o,^^^-^^^2-^i^^iJXO _ ^324) 



{H^-H[)xo = XiHoXo (B.22) 
{H,-H[-\,Ho)xi = -iH2-H',-\2Ho-\iH[)xo. (B.23) 

Note that the zeroth order is fully degenerate and Xo is determined only 
in first order by again a generalized eigenvalue problem. This is similar to 
ordinary highly degenerate perturbation theory, where a large diagonalization 
cannot be avoided. Combining the first two orders gives also 

xl{H2 -H!,-\ ,H[)xo 
XoHqXo 

Two generalized eigenvectors x, x of ()B.1|) belonging to different general- 
ized eigenvalues A 7^ A are orthogonal in the sense 

= x^Hx = X^H'x. (B.25) 

In particular, eigenvectors to be constructed now have to be orthogonal to 
H or H' times the two previous ones. Of course, this has to come out 
automatically, but it is interesting to verify as a consistency check. Using 
the coincidence with eigenvectors 0,0' up to terms of 0{g^), the space, to 
which Xo must be orthogonal, may be spanned by the leading orders of Hcf)' 
and H'(f). A little calculation in ordinary perturbation theory gives 

H<f)'-H'<p = 2g{H,-H[)<Po + 0{g') (B.26) 
Hcl>' + H'cl> = g\H,-H[)H,\H,~H[)(l)o + 0{g^). (B.27) 

The orthogonality of solutions xo to these directions follows indeed by first 
projecting both sides of ()B.22|) on 0q and then on (j)\{Hi — H[)Hq^ and using 
flR^ again. 
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In our real problem we have to deal with the additional spin degree of 
freedom. Hence there are two zero-modes (f)oa,<y = 1,2 of h^. These are 
chosen orthonormal and such that the 2x2 matrix 

0L(^2 - hi{hoy^hi)(j)of3 = K2a 5af3 (B.28) 

is diagonal. The expansions built upon them generate two generalized eigen- 
values which are 0{g^) and two 0((?~^). 



B.2 Evaluation for given gaugefields 

Here we evaluate ^2, A2, A'„2 for gaugefields 

A^{x) = yJ2 ^m(p) exp(ip ■ x) 



(B.29) 



and an analogous expression for A'. According to ()2.15|) we have (for gauge 
parameter = 0) 

A,{p) = Y.^,Al-e"'^^)Mp) (B.30) 



with 



{B.31] 



leading to a real A^. 

We consider the expansion of h = Dw'j5 now in the ON-basis of free field 
states (l,0)e*^^'/L and (0, l)e*^^'/L, i.e. as kernels in momentum space and 

o 

2x2 matrices in spin. With = 2 sin(p^/2),P^= sin(p^) we find 



ho{p,q) 
hi{p,q) 
hip^q) 



PI 



;b.32) 



^ Y: [il, + l)e-^^'' + (7;. - l)e^^''] ^^A.ip - q) (B.33) 

At 



+ l)e-*^'- - (7^ - l)e*^- 75^%(P - q) (B.34) 



with A'^^{p — q) = J2r Af_i{p — r)A^{r — q). Combinations relevant for ()B.28|) 



are 



/i2(0,0) 



75 



2L2 ^ 



T.p"\Mp)\' 



hiip,o) = 7rFY.^i^^f^(p)Mp) 



(B.35) 
(B.36) 
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with 



(zo, zi) = -2te-'^P°+P'^/^ (pocos(pi/2) , pi cos(po/2) 

With some calculation we are then able to show 

h{0,p){ho{p,p))-'h{p,0) = 
-1 



4L2 



75 



(B.37) 



(B.38) 

\Up)? 



Upon p-summation the part odd in p does not contribute and from ()B.28j) 
we get 

1 



^22 — —1^21 



2L2 



E 



1 + 



"2 1 "2 "2 

P - WoPi 

p +Hp'r 



p \Mp)\^ 



(B.39) 



Inspecting fl2.18|l we find on average in the quenched ensemble 

1 



(IMpW) 



G 



A numerical asymptotic expansion gives 



K22 )g = — In L - 0.00639327 + 0.13303 L^^ + 0{L-^) 
2tx 



(B.40) 



:B.41) 



Also the mean of the square can be worked out 

( {ti2af )g = { H2a )g + 0.00773389 + 0.03112 L'^ + 0{L-^). (B.42) 

The denominator in ()B.17|) is given by elements of the 2x2 matrix 

J2ih[ip, 0) - hip, 0)yih[ip, 0) - hip, 0)) = (B.43) 
p 
1 



j2T.ip' - 2^1pI) \€ip)-Mp)\' 



and produces 



( -/ii)0oar )g = -hiL - 0.06162368 + 0.26607 L-2 + 0(L-^). (B.44) 

TT 
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